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ABSTRACT 

/~The  integro-differential  equation  describing  the  diffusion 
of  neutrons  in  a  multiplying  medium  is  studied  under  the  assumption 
of  plane  geometry  and  linearly  anisotropic  scattering,  hut  without 
the  usual  limitation  to  continuous  slowing  down  and  weak  absorption. 
Under  the  condition  that  the  integrated  flux  vanish  at  the  boundary, 
it  is  shown  that  a  separable  solution  exists.  The  relation  between 
core  size  and  the  number  of  neutrons  required  per  fission  for  a 
steady  state  is  determined.  The  results  are  specialized  to  the 
case  of  isotropic  scattering  in  the  laboratory  system,  and  to  the 
case  of  a  medium  containing  hydrogen  and  an  element  which  scatters 
without  energy  loss.  /  a  U-  A (U  J 
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I.  INTRODUCTION 

In  designing  a  nuclear  reactor,  one  of  the  major  considerations 
is  the  amount  of  fissionable  material  required  to  sustain  a  chain  re¬ 
action,  If  neutrons  are  moderated  by  elements  of  intermediate  atomic 
weight  for  which  inelastic  scattering  is  negligible,  and  if  little  ab¬ 
sorption  is  present,  the  problem  can  be  treated  by  age  -  diffusion 
theory.  Two  cases  of  interest,  however,  do  not  fall  under  this  formal¬ 
ism:  these  are  moderation  by  a  heavy  element,  where  energy  losses  are 
caused  primarily  by  inelastic  scattering,  and  moderation  by  hydrogen, 
in  which  cross-sections  change  radically  in  a  single  collision  interval. 

The  purpose  of  this  report  is  to  study  these  two  cases  in  terms 
of  the  integro-differential  equation  underlying  diffusion  theory.  The 
discussion  will  be  limited  by  the  following  assumptions:  First,  the 
reactor  has  plane  symmetry,  so  that  the  neutron  distribution  depends 
only  upon  a  single  space  coordinate.  Second,  the  probability  that  a 
neutron  will  be  scattered  through  a  given  angle  can  be  represented  by 
a  constant  plus  a  term  proportional  to  the  cosine  of  the  angle.  Finally, 
the  neutron  density  at  the  boundaries  is  zero  for  all  energies. 
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II.  ASYMPTOTIC  SOLUTION  OF  THE  BOLTZMANN  EQUATION 

We  asstime  the  reactor  consists  of  a  slab  bounded  by  vacuum,  ex¬ 
tending  along  the  x-axis  from  -a  to  +a  and  indefinitely  in  the  y-z  plane. 
Its  time  -  independent  state  is  then  described  by  the  function  N(x,v,^x) 
dx  dv  dju,  defined  as  the  number  of  neutrons  in  a  volume  element  'with 
unit  cross-section  and  thickness  dx,  having  scalar  velocities  between  v 
and  v  +  dv,  and  with  velocity  vectors  making  an  angle  with  the  positive 
x-axis  whose  cosine  lies  between  jp.  and^x  +  dju.  It  is  slightly  more  con¬ 
venient,  however,  to  take  as  a  distribution  function 
F  -  N  U<VVAA)  V  , 

Two  quantities  needed  later  can  be  defined  in  terms  of  F;  these  are  the 


(2.1) 


and  the  current 


F  i  (*1*0  ~ 


(2.2) 


The  form  of  the  continuity  equation  we  shall  use  requires  that  the 
differential  scattering  cross-section  be  linear  in  the  scattering  angle. 
The  differential  scattering  cross-section  0^(v,v',^i0)  dv  d^iQ  is  defined 
as  the  macroscopic  cross-section  for  scattering  from  the  scalar  velocity 
v  into  the  interval  from  v  to  v  +  dv,  through  an  angle  (measured  in 
the  laboratory  system  of  coordinates)  whose  cosine  lies  between  and 
Po  +  ^o*  &eneral  ^(v,v/,^xq)  can  be  expanded  in  an  infinite  series 
of  Legendre  polynominals  in  We  will  make  the  assumption  that  only 

the  first  two  terms  are  different  from  zero,  and  write 

cr^  (V.V'.yUo)  =  <rsw>  ■+ 


Using  the  orthogoality  of  the  Legendre  polynominals,  and  integrating 
this  expression  with  respect  to  it  follows  immediately  that 


^  (V,  v") 


=  Os  Cv.V'.^c) 


and  is  the  cross-section  for  scattering  from  v  to  a  unit  interval 

about  v  through  any  angle,  while 

_  f  VV**  CT^  (v.  V',  A,) 

ueiv)  =  - — 

'  jV^CV.V',/6) 

and  is  the  average  cosine  of  the  angle  through  which  a  neutron  is 
scattered  in  a  single  collision. 

Under  the  assumption  that  neutron-neutron  interactions  are  negli¬ 
gible  the  distribution  function  F(x,v,yu)  satisfies  the  following 
equation  of  continuity  due  to  Boltzmann:  ^ 

aF(*/VO  ^  f  (%,*,/».}  -  {  %)  Xm  \  A/*'  nx,vj/0 

'  Jo 

4  i  f^v'Ts  (V,v'>  fv'  F  tx<v///^1)  f  |  ^  feW^tv'ja^v.v') 

2*  JV  v  *  i 

where  cr(v)  and  c^(v)  are  the  macroscopic  total  and  fission  cross- 

sections  respectively  for  neutrons  of  velocity  v;  }C(v)dv  is 
fraction  of  fission  neutrons  born  with  velocities  between  v  and 
v  +  dv,  vQ  is  some  finite  velocity  above  which  no  fission  neutrons 
are  born,  and^is  the  number  of  neutrons  produced  per  fission. 

It  is  assumed  that  both  i)  and  %  are  independent  of  the  energy  of  the 
neutrons  causing  fission,  and  that  the  angular  distribution  of 
fission  neutrons  is  isotropic  in  the  laboratory  system. 

A  detailed  derivation  of  equation  (2. 3)  is  given  by  Weinberg"*" 
whose  notation  is  largely  followed  here.  The  two  terms  on  the  left 


XA.  M.  Weinberg,  Pile  Neutron  Physics,  M-3336,  equation  1-28. 
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represent  the  loss  of  neutrons  from  a  unit  volume  due  to  leakage 
through  the  boundaries  and  to  absorption  and  scattering  within  it. 

The  first  term  on  the  right  corresponds  to  neutrons  produced  by- 
fission  and  the  second  to  neutrons  produced  by  scattering  down  from 
a  higher  energy.  The  last  term  represents  a  distortion  in  the  angu¬ 
lar  distribution  of  these  scattered  neutrons  due  to  the  anisotropic 
scattering  in  the  laboratory  system;  if  -  0,  this  term  is  absent. 

The  cross-sections  e r  ,  and  cr ,  through  which  the  dependence 

of  F  on  the  material  properties  of  the  reactor  enters,  have  been 
written  as  dependent  only  on  velocity.  However,  spatial  dependence 
can  be  taken  into  account  provided  it  is  of  the  form 

crc*,v>  ~  pcx)  c r  w) 

where  cr  (x,v)  refers  to  any  of  the  cross-sections  and  p(x)  is  some 

2 

fixed  function  of  X.  If  the  new  coordinate 


is  introduced,  together  with  the  function 

F'(  X'OO  /  v, yA.)  ~  F 

p(x)  can  be  cancelled,  giving  an  equation  in  F*  formally  identical 
with  (2.3) •  It  is  therefore  sufficient  to  consider  (2.3)  as  it 
stands . 


This  transformation  is  a  special  case  of  the  optical  distance 
used  by  B.  Davison,  Transport  Theory  of  Neutrons,  LT-18,  Chalk 
River,  Ontario,  Rational  Research  Council  of  Canada,  January  19^7> 
p.  10. 
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The  problem  we  consider  is  the  following:  knowing  the  values 
of  the  cross-sections  and  of  for  a  particular  composition,  we  re¬ 
quire  the  critical  size  of  the  reactor,  that  is,  the  size  for  which 
a  steady  state  will  he  maintained.  Physical  considerations  indicate 
that  a  solution  does  not  always  exist,  since  enough  parasitic  ab¬ 
sorption  may  be  present  to  prevent  a  self-sustaining  chain  reaction 
under  any  conditions.  It  is  therefore  more  convenient  to  consider 
the  size  as  given,  and  to  determine  the  value  of  t>  required  for 
criticality.  We  accordingly  take  the  boundaries  at  x  =  +  a  as 
given  and  look  for  a  solution  of  (2-3)  for  F(x,v,ju)  together  with 
the  corresponding  value  of  ,  subject  to  the  conditions 
pfx, 3  F  (-  X  .  v,  -/•*-) 

>  o  icr  -a  <  x  <  a 
/=0ca.v')  =  F„  =  O  ior  all  v 

The  first  condition  is  due  to  the  symmetry  about  x  =  0,  and 
the  second  follows  from  the  physical  requirement  of  a  positive  neu¬ 
tron  density  within  the  reactor.  The  last  condition  is  only  an 
approximation  to  the  physical  system  and  is  introduced  (as  is  usually 
done  when  considering  energy  dependence)  to  simplify  the  problem. 

The  proper  boundary  conditions  for  a  reactor  bounded  by  vacuum  are 


F 

)  =  o 

■f  o  v 

-  |  4L  <  O 

F 

(-a  .v' 

y*)  E  O 

{  CV 

o  <  /*>  <  1 

corresponding  to  no  incoming  neutrons.  Under  these  conditions,  the 
space  distribution  of  the  flux  has  been  determined  for  the  case  of 
monoenergetic  neutrons  (ref.  3),  and  is  found  to  approach  asymptoti¬ 
cally  the  distribution  obtained  using  the  approximate  condition  (2.6), 


(2A) 

(2.5) 

(2.6) 
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the  difference  between  the  two  distributions  decreasing  away  from 
the  boundary  like  an  exponential  with  a  relaxation  length  of  a 
mean  free  path.  It  seems  reasonable  to  expect  a  similar  behavior 
for  the  energy  dependent  case;  the  error  in  critical  size  caused 
by  using  the  approximate  boundary  condition  should  then  be  small 
if  the  reactor  thickness  is  large  compared  to  a  mean  free  path, 
as  it  usually  is. 

In  the  remainder  of  this  section  ve  will  show,  following  a 
procedure  used  by  Weinberg^  for  the  case  of  monoenergetic  neutrons, 
that  solutions  of  the  energy  dependent  Boltzmann  equation  exist 
having  a  sinusoidal  dependence  on  x.  Since  the  equation  is  linear, 
it  is  convenient  to  work  with  the  complex  function 

&  *  e‘  j  tv,/*.)  ■  (2.7) 


whose  real  and  imaginary  parts  will  both  satisfy  (2.3)  if  F  itself 
does.  We  make  the  tentative  assumption  that  a  solution  of  the  form 
(2*7)  exists  and  try  to  determine  the  function  f  (v,yu)  and  the  para¬ 
meter  K.  Substituting  this  in  place  of  F  in  (2-3)  and-  cancelling  the 


exponential  leads  to  i 

(o-rcyo-fovo  =  f  iv.^cvo^'fcvy.) 


(2.8) 


jl  |  +  |/«.J  av'^voo^v.v'J 


(2.9) 
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i 

>  S  l  f 


2.10) 


CvO  a  V 


Xcv)  I  ^v'0^(v'>  f0(v')  +  J  Av’crs  cv.v') /oCV) 


2.11) 


aw  tv 


3  f  ^V'yA.tw')  0*s  W,V») 


2.12) 


equation  (2.8)  can  "be  written 

f  ,  I  /-H tv) 


5  <-v/o 


X  0>)  +  JK/W. 


2.13) 


Substituting  this  expression  for  f (v,ju)  in  (2.9)  and  (2.10) 
using  the  integrals 


cK/^  . 

i^svT 1 

'  04 

-  —  5  n 

-  a  Oi 

-  A/3 

=  /. . 

i*vsT*c*  \ 

1-/3- 

l  4-LOyA. 

lot  \ 

t 

o< 

=  A  /, 

_  \ 

r  2 

i-  r> 

\  +  loy* 

C*.L  k  1 

04.  / 

CK  x 

(2.1k) 


K  0“  and  p  =  /o4 

,  we  have 

ie(V)  -  £  -  1 

06  / 

(2.15) 

■f  Cv)  -  —  (  2  AA. 

J  <r  v  ot 

-b  1  ~  —  ^vw  ) 

«>  ' 

(2.16) 

it  being  understood  that  m,  n,  06,  p  and  cr  are  velocity  dependent 
and  that  tan^oc  is  restricted  to  its  principal  value.  When  these 
expressions  for  fQ  and  { ,  are  substituted  in  (2.11)  and  (2.12),  the 
following  simultaneous  integral  equations  for  m  and  n  are  obtained 


/yyt(v)  - 


d-v'  sn(v‘)  -t  f  °cL v1  KU\v,v')  *\(v')  (2.17) 
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MtV)  5  2  (v 


)  _  ( CUV  KflV,  V')  f**  tv<)  +  ^v'Jc^.v') /n(y') 

Jv  Jv  (2. 


(2.18) 


where 


_  .  /  —  i^s C v 0  0%  (V,VO 

*  -  v**  «•■>  -^r 

Kmiv,v')  =  3A.  (V)  1  ~ h t'".)  <U  o.i'O 

'  a(v')  —  - 


<Tiv> ) 


k<3)(*,v o  =  ft q~s  (V/ vo 
cy(v")  crevo 

/C<4)  ry.vO  _  ^(v#)  <r$  (w,  v>»; 

CTO') 


Q.(v)  =  vXh)  S-H  2  f/^CV')  <n(V')  +  iz£±¥2  'VvUV") 

J°  V(V)  L  cx'(V') 


Xtv) 


(2.19) 


c  being  an  undetermined  constant. 

In  order  to  show  that  (2.17)  and  (2.18)  actually  define  the 

1*. 

functions  n  and  m,  we  require  the  following  theorem: 

The  integral  equation 


At  tv) 


f  V'o 

_  frfv)  +  J  dv>  K(v,V)  K(V') 


(2.20) 


where  |K(v,v/)|  <.  M  for  all  v,  v*,  has  exactly  one  solution  which  is 


finite  everywhere.  This  solution  is  given  by 

rv° 

At-tvO  s  g  (v)  +  J  Av>  N(V,V)  2  O’) 


(2.21) 


where 


N  tv,  v1 ) 


K^<W) 


is  the  Neumann  series  in  the  iterated  kernels,  the  latter  being  de- 


4  t 

V.  Volt err a,  Lecons  sur  les  Equations  Intfegrales,  Paris,  Gauthier 
Villars,  1913,  p.  50  -  52- 


fined  by 


■K,  Cv,v')  =  K  Cv.v') 

V  1 

X  (V, V)  =  K  (VJW)  (2.22) 


The  solution  of  (2.17)  and  (2.18)  can  now  be  obtained  by 
a  procedure  similar  to  that  used  for  simultaneous  algebraic 
equations.  Using  the  theorem  just  proved  and  considering  n(v) 
as  given,  equation  (2.17)  can  solved  for  xn(v)  in  terms  of 
n(v): 


r  ° 

•  jy 


r*(V) 


,v i  ,*>  . 

+.  7v',v")  l 

Ju  JY 


(2.23) 


■(1) 


where  (v,v')  is  the  Neumann  series  corresponding  to  K''* '(V^V*) 

Substituting  this  expression  for  m(v)  in  (2.18)  yields 


an  equation  of  the  form 


/n(v) 


g.(v)  +  J  *L‘"  ^Cvl) 


{2.2k) 


containing  n(v)  only,  where  g(v)  is  the  quantity  defined  by  (2.19) 
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Equation  (2.210  has  the  solution 


mLV)  -  ^  + 


ra 

l  N(V,  V)  ^  (V) 


(2.26) 


where  N(v,t?)  is  the  Neumann  series  in  the  kernel  (2.25)  and  m(v) 
can  now  be  found  from  (2.23).  m(v)  and  n(v)  as  defined  by 

(2.23)  and  (2.26)  are  a  solution  of  the  system  (2.17)  and  (2.18) 
is  evident  since  all  operations  used  in  deriving  (2.23)  and  (2.26) 
are  reversible.  The  pair  of  functions  m(v)  and  n(v)  determined  by 
(2.17)  and  (2.18)  are  unique  for  a  given  g(v)  since  a  different 
pair  would  lead  to  the  same  equations  (2.23)  and  (2.26) . 

In  order  to  determine  V,  the  values  of  n(v)  and  m(v)  so 
obtained  are  substituted  in  the  equation 


(V) 


-h 


aUv) 


rf'VVV  (.V 


c 


(2.27) 


which  defines  c  (cf.  (2.19)).  It  is  clear  from  the  form  of  (2.23) 
and  (2.26)  that  n(v)  and  m(v)  will  both  contain  c  as  a  factor;  it 
can  therefore  be  cancelled  from  (2.27)  leaving  an  equation  in  which 
all  quantities  are  known  except  v* .  The  presence  of  the  arbitrary 
constant  c  is  caused  by  not  specifying  the  absolute  value  of  the 
flux,  and  corresponds  to  the  fact  that  a  critical  reactor  will 
maintain  a  steady  state  at  any  power  level.  Since  the  value  of  c 
has  no  effect  on  the  functional  form  of  m(v)  and  n(v)  or  the  value 
of  V,  it  will  be  taken  as  unity,  so  that  m  and  n  are  both  real. 

The  complex  solution  of  the  Boltzmann  equation  can  now  be 
written,  using  (2.7)  and  (2.13),  as 


I  Alt*)  -C/A/v*.(W 
g.  =  acTft-T  i  +  i  o(v\/a.  ~ 

I  (/vi  ,  «,«*'*•  Av.  ),  Ctf-o  K  X  4-  (o<  M  +  '’Vn.  )  K  X 

acr  i  +  o<.^  /*.•«• 

L  C/n  -  otyu  tAW,  )  /a^U.  k  x  — /*■  (qi<ia  +  KX 

Z<r  i  v  c*-*-^*- 


(2.28) 


That  this  function,  with  m(v)  and  n(v)  defined  hy  the  integral 
equations  (2.17)  an*3-  (2.l8),  is  a  solution,  can  be  seen  by  sub' 
stituting  it  in  place  of  F  in  (2.3) •  On  cancelling  the  ex¬ 


ponential 


and  equating  coefficients  of  like  powers  of  p, 


the  identities  (2.11)  and  (2.12)  result.  Both  the  real  and 
imaginary  parts  of  F  are  solutions  of  (2.3) >  however,  only  the 
real  part  satisfies  the  symmetry  condition  (2.k).  Using  the  de¬ 
finitions  (2.1),  (2.2),  (2.9)  and  (2.10),  the  integrals  (2.1k) 

and  retaining  only  the  real  parts, 

„  .  N  _  ^  c«KK  +/*  (o + 'Vw  )  Au*- K  * 

“  *«-  (2.29, 

Ft,  (X>V)  -  (fi  av  -  sxvx.')  Qjo  Ac*  -  io^  OrofcX  (2*30) 

F,  t*.*)  =i  (■*§£«■+  !g~)^KX  (2-31> 

Conditions  (2.5)  and  (2.6)  can  be  met  by  choosing 


(2.32) 


thereby  determining  K  in  terms  of  the  reactor  dimensions. 

Most  of  the  previously  defined  quantities  can  now  be 
interpreted.  From  (2.30)  and  (2-31) >  fo  ij  are  the  enerSy 
dependent  factors  of  the  flux  and  current  respectively.  From 
(2.8),  n(v)dv  is,  neglecting  the  space  factor,  the  number  of 
neutrons  entering  the  velocity  interval  dv  in  unit  time  by  being 
bom  there  as  fission  neutrons  or  scattering  into  the  interval 
from  a  higher  velocity.  Since  a  steady  state  exists,  it  is  also 
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the  rate  at  which  neutrons  are  lost  by  absorption,  scattering  and 
leakage.  Comparing  (2.12)  with  (2.3),  it  can  be  seen  that  m(v) 
is  essentially  the  distortion  term  due  to  anisotropic  scattering; 
its  effect  ((2.25)  and  (2.26))  is  to  decrease  the  flux  and  increase 
the  outward  current.  The  quantity  g(v)  (2.19)  is  the  fission  source 
term.  Finally,  ot  =  K  /  (T  is,  using  (2.32)  and  neglecting  the  factor 
jt/2,  the  ratio  of  the  total  mean  free  path  to  the  half -width  of  the 
reactor,  a  quantity  assumed  to  be  snail. 

In  satisfying  the  condition  that  the  flux  vanish  on  the 
boundary,  it  would  be  sufficient  to  take  the  parameter  k  equal  to 
any  odd  multiple  of  jt/2a,  so  that  the  solution  does  not  appear  to 
be  determined  completely.  For  the  case  of  isotropic  scattering  in 
the  laboratory  system,  however,  we  shall  show  that  only  (2.32) 
leads  to  a  solution  of  the  problem.  By  properly  choosing  the  com- 
plex  functions  j  (v,^i),  any  solution  of  (2.3),  (2.k),  (2-5)  and 
(2.6)  can  be  written  as  the  real  part  of 


~  M)  iK  «* 

si-o 


(2.33) 


where  (  3  +  i) 


(2.3*0 


Substituting  (2*33)  in  (2.3)  and  equating  coefficients  of  & 
leads  to  an  equation  of  the  form  of  (2.8)  for  each  f ^  .  We 
suppose  is  a  particular  solution  of  the  problem  with  t> 

the  corresponding  number  of  neutrons  per  fission  required  for 

t  ,  M 

criticality,  and  consider  any  two  distinct  modes,  say  £  and  j  . 
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Following  the  previous  procedure  with  -  0  leads  to  the 
equations,  corresponding  to  (2.18)  and  (2*32), 


r  X  <*)  +  \  tftV  K  (V,  V*')  /H'(V') 

J  y 

2  _  f  dv  Sk. /a'cvj/n'fv) 

P  “  Jo  OW)  1 


(2.35) 


„>/(v)  =  X(v)  + 

fVo 

J/v' 

r^\J 

j-  ^  [  dv  q?tvj 

(i'Vv) 

Av"  (vO 

J»  <r  iv) 

/w"lv') 


(2.36) 


where 


cr^vO 


The  constants  c*  and  c**  have  both  been  taken  as  unity  since  only 
is  of  interest.  Suppose  K,;  >  K*  •  Then  o<.  “  =  <  /CT  >  ^ 

-I  .  <1  it 

and  B"  «  _  <  pf  for  all  v.  The  Neumann  series  in  K  (V,V0  is 

ot  “ 

accordingly  less  than  the  one  in  K  iV|V')  and  since 


^’iv)  •  "Xcv)  +■  f  *v>  Xiv>) 

'  TCM  +  V(V/)  X<*"> 


hy  (2.21),  it  follows  that  n"(v)<  n'  (v)  for  all  v. 

From  (2.35)  an13-  (2.3^), 

f\v  °^(v)  r  fi> 'M  S*'(v)  -/*"(*)  S»"(V)1  ^l-±  >  © 

4  o*W  1-  J 

In  view  of  this  contradiction,  the  initial  assumption  that 
the  solution  contained  more  than  one  mode  cannot  be  true;  if 
several  solutions  of  the  problem  exist,  they  each  consist  of 
only  a  single  mode.  But  the  only  such  solution  which  is  never 
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negative  in  the  reactor  is  the  fundamental  ( j  =  o  in  (2.3^-))>  which 
therefore  constitutes  the  unique  solution  of  the  problem. 
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III.  ISOTROPIC  SCATTERING 


In  this  section  some  of  the  previous  results  will  he 
specialized  to  the  case  in  which  scattering  is  isotropic  in  the 
laboratory  system,  and  will  apply  approximately  to  a  reactor 
consisting  of  heavy  elements  where  moderation  is  due  almost  en¬ 
tirely  to  inelastic  scattering. 

Since  ^  and  hence  m(v)  are  zero,  the  flux  and  current  are 


given  by 


F  C* 


,v)  =  Jo  (V)  ^  K  X  * 


ft  tv) 


F,  cv)  KX 


J  "  *v(y  )  A**-  K  X 

crtv)ottv) 


(3*2) 


•Where  n(v)  is  the  solution  of  the  single  equation 


/ 

=  X‘»  *  J  dv  -W)  *  <v') 


/ntv)  - 


end  V  is  then  determined  from 


(3-3) 


?>  (  av  /uv)  ~  i 

J»  <r(v) 


(3-4) 


corresponding  to  equations  (2*30) >  (2*31) >  (2.18)  and  (2.27) 


respectively. 

It  is  convenient  to  introduce  the  following  notation: 

The  probability  that  a  neutron  which  has  left  velocity  v  by 
any  means  (scattering,  leakage  or  absorption)  will  be  scattered 
to  a  velocity  between  v  and  v  +  dv,  will  be  designated  by 
a( v,  v#)  dv.  Since  this  is  the  ratio  of  the  rate  at  which  neu¬ 
trons  leave  v  by  scattering  into  the  interval  dv,  to  the  rate 


at  which  they  leave  by  any  means, 
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v  Os(v.V)  £tv')  Gn>  k*  _ Os(v,V') 

dW,v)  =  -  =  fiOfO  :■ 

/VUV  )  Oro  K  X  u  tV) 


(3-5) 


The  last  equality  follows  from  (3.1),  which  can  he  written 


Similarly,  4 [v)  is  defined  as  the  probability  that  a  neutron 
which  has  left  velocity  v  has 'done  so  by  causing  fission,  and 


is  given  by 


(v)  r  ^CV)  ^  KX  -  /MV)  CV) 

oacv)  Oro  K  %.  (T  <V) 


(3.6) 


Using  these  definitions,  (3*3)  and  (3*^)  can  be  written 


/ruv)  -  X  M  +  1  <Av'alv,v*>  /k(W) 


(3-7) 


V  fo**)  -  I 

J  Ct 


(3-8) 


The  Neumann  series  solution  of  (3*7)  can  be  written 


(v)  *  X  ^  iv) 


(3-9) 


provided 


/n#  cv)  -  X 


I 

ty)  e  <&V,a.LVlV')  M-j  (V) 
\+i  • 

v 


(3.10) 


as  can  be  seen  by  comparing  (3*9)  and  (3*10)  with  (2*21)  and  (2.22) 
From  (3*10)  the  quantity  n.(v)dv  can  be  interpreted  as  the  number 
of  neutrons  entering  the  interval  dv  which  have  been  scattered  j 
times  since  birth  as  fission  neutrons.  It  is  evident  that  this 
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method  of  solution  is  not  well  suited  to  numerical  computation 
unless  such  heavy  absorption  or  leakage  is  present  that  a  neu¬ 
tron  makes  only  a  few  collisions  before  being  lost  from  the  re¬ 
actor  .  An  alternative  is  to  start  with  a  function  more  nearly 
representing  the  final  distribution  than  the  fission  spectrum 
and  proceed  by  iteration.  However,  the  kernel  does  not  have  a 

i 

simple  form  and  the  integrations  must  be  carried  out  numerically. 
It  is  therefore  desirable  to  replace  (3*7)  by  a  se^  simulta¬ 
neous  algebraic  equations  which,  because  of  their  forms,  can  be 
solved  immediately.  If  the  velocity  range  is  divided  into  N 
intervals  numbered  in  order  of  decreasing  velocity,  (3*7)  and 


(3.8)  can  be  written 


;  -  X;  +  2  aij 
0-  1 


(3.11) 


(3.12) 


where  n^  is  the  number  of  neutrons  leaving  the  i  interval  in 
unit  time,  ^  is  the  probability  that  a  neutron  leaving  the  i 
interval  has  done  so  by  causing  a  fission, 

^  is  the  fraction  of  fission  neutrons  born  in  the  i^11 
interval,  and 

a,  .  is  the  probability  that  a  neutron  leaving  interval  j 

J.  J 

has  been  scattered  to  interval  i.  Since  a  neutron  which  has  been 
scattered  into  the  same  interval  is  considered  not  to  have  left 
that  interval,  a^  =  0  and  the  upper  limit  on  the  summation  in 
(3.11)  is  i-1.  The  n;-  's  can  be  found  recursively,  starting  with 
the  highest  velocity: 
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*-*l  +  /H' 

=  Xx  +  a3t  +  a3l  ^ 


(3-13) 


after  which  Vis  found  from  (3*12) 

Instead  of  using  the  continuity  equation  (3*7)>  one  can 
work  with  the  following  adjoint  equation: 


=  +  [ 


(LV>  CL(Vf,V  )  fCY1) 


(3-1^) 


which  defines  the  function  p(v)  in  terms  of  the  known  quantities 
(v)  and  aCv^v).  In  order  to  find  a  condition  similar  to  (3-8) 
which  will  determine  V  in  terms  of  p(v),  multiply  (3*7)  ty  p(v) 
and  (3.14)  by  n(v)  and  integrate  with  respect  to  v  from  o  to 

j  °aiY  p(v)  s  J f>cv)  X(y)  4  ^  <A v  piv)  J ctY'eUv,i',)A'(v9 

J  dv  ?  tY)  r*.tv)  -  J  +  j  olv'  AWi/)  PCv') 

The  two  double  integrals  are  equal,  as  can  be  seen  by  reversing 
the  order  of  integration  in  either  one  and  interchanging  v  and 
v ' .  Subtracting  one  equation  from  the  other,  and  making  use  of 


(3*8),  we  have 


v 

|  JLv  ptY) 


XCv)  =  / 


(3-15) 


as  the  required  condition, 


One  can  therefore  interpret  p(v)  as  the  probability, 


averaged  over  the  reactor,  that  a  neutron  with  velocity  v  will 
eventually  result  in  a  fission.  Equation  (3*1*0  is  the  state¬ 
ment  that  this  probability  is  the  sum  of  two  parts;  the  pro- 
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bability  of  causing  fission  immediately  on  leaving  v;  and  the  pro¬ 
bability  of  first  being  scattered  to  some  lower  velocity  v  7  , 
times  the  probability  that  a  neutron  at  v'  will  eventually  cause 
fission,  summed  over  all  v7.  Under  this  interpretation  both  (3°1*0 
and  (3-15),  which  requires  the  neutron  distribution  to  reproduce 
itself,  could  have  been  written  immediately,  but  the  explicit  form 
of  a(v,v')  would  not  have  been  apparent. 

For  numerical  solution,  (3*1*0  and  (3*15)  can  be  replaced  by 


and 


(3-17) 


The  solution  of  (3*l6)  is 


1°N  * 

Pm-1  =  £a/-i 


+  A NN-I 


(3-18) 


Formulating  the  problem  in  this  way  has  some  advantage  over  the 
continuity  equation,  since  the  fission  spectrum  enters  only  after 
the  p^*s  are  determined.  This  permits  investigation  of  the  effect 
on  criticality  of  various  assumptions  for  the  fission  spectrum 
without  repeating  the  major  part  of  the  calculations. 

The  adjoint  equation  is  closely  related  to  the  method  of 
successive  collisions  of  De  Marcus  and  Nelson  (ref.  2)  which  they 
have  used  to  study  the  transmission  of  particles  through  matter. 
Their  method  consists  of  describing  the  neutron  distribution  by 


a  state  vector  whose  component  is  the  number  of  neutrons  in 
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the  i  th  cell  in  phase  space,  each  cell  corresponding  to  a 
particular  specification  of  position,  velocity  and  direction, 
of  which  there  are  a  finite  number.  In  addition  to  these  'free 
states1,  ’trap  states’  are  included,  corresponding  to  processes 
which  remove  neutrons  from  the  free  states,  and  having  the  pro¬ 
perty  that  a  neutron  cannot  leave  them.  A  matrix  operator  A  is 
introduced  whose  element  A- i  is  the  probability  that  a  neutron 
which  has  left  state  j  will  go  to  state  i  (which  may  he  a  trap 
state) .  After  k  transitions,  the  neutron  distribution  is 

s-  /\  ^  (p 

The  ultimate  distribution  is  given  by 

=  A*  V  (3-19) 

where 

i*  = 

If  this  limit  exists,  the  equation 

A*  A  •  A"  (3'20> 

holds.  Since  A  is  known,  A*  can  be  determined  from  (3.20),  and 
then  used  in  (3*19)  to  find  the  ultimate  distribution  of  neutrons 
among  the  trap  states. 

Although  this  method  was  set  up  for  a  non-reproducing  medium, 
it  can  be  used  to  find  the  critical  size  of  a  reactor  by  taking  as 
the  trap  states  capture  leading  to  fission,  using  the  fission  spect¬ 
rum  for  the  initial  distribution,  and  considering  only  one  neutron 
generation.  A  generation  is  taken  as  the  interval  between  birth  as 


a  fission  neutron  and  death  by  absorption  or  leakage.  Under  the 
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conditions  we  have  assumed,  the  application  of  this  method  is  greatly- 
simplified,  since  the  space  distribution  is  known  in  advance.  This 
reduces  the  number  of  unknowns  and  makes  it  possible  to  solve  the 
resulting  simultaneous  equations  recursively. 

We  take  N  free  states,  each  state  corresponding  to  a  particular 
velocity  range  and  numbered  from  1  to  N  in  order  of  decreasing 
velocity.  In  addition,  we  take  the  N+l  st  state  as  a  trap  state, 
to  which  neutrons  causing  fission  are  sent.  Making  use  of  its 
definition,  the  matrix  A  is  then  of  the  form 


( 


o  o  o  -  *  *  o  c 

a2l  °  O  '  '  ‘  0  O 

a3*  a31  O'1'  O  0 

*  *  •  • 

aN9  '  •  *  O  o 

fi  n  •  •  *  i**  I 


(3.21) 


The  sub-matrix  consisting  of  the  a.  . ’s  is  triangular  since  neutrons 
are  only  scattered  downward  in  energy.  The  last  column  has  the  form 
shown,  since  neutrons  in  the  trap  state  cannot  leave  it,  and  the  form 
of  the  last  row  follows  from  the  previous  definition  of  the  f±'s. 


It  should  be  noted  that  the  sums  of  the  elements  in  each  column  are 
generally  less  than  one.  If  additional  trap  states  were  included, 
corresponding  to  other  ways  by  which  neutrons  leave  the  reactor,  it 
would  be  possible  to  account  for  all  the  neutrons  involved.  For 
our  purpose,  however,  only  the  disposition  of  those  neutrons  causing 
fission  is  of  interest. 

Assume  a  finite  probability  of  loss  in  each  transition  from  a 
free  state,  and  consider  the  probability  of  passing  from  an  initial 
free  state  j  to  some  final  state  after  an  infinite  number  of  transi- 


27 


tions-  This  is  clearly  zero  if  the  final  state  is  a  free  state, 
and  pj  (cf.  the  definition  on  p.23)  the  final  state  is  the 
fission  trap  state.  The  iterated  matrix  is  therefore 


(3.22) 


Substituting  (3«2l)  and  (3*22)  in  (3*20)  and  multiplying  out 
leads  to 


+ 


(3.23) 


If  we  take  as  the  initial,  distribution 


the  ultimate  disposition  of  the  neutrons  of  a  particular  genera¬ 
tion  will  be  given  by  the  N+l  st  component  of  the  column  vector 
A*X  which  is 

5  t1:  X  ; 

In  order  for  this  generation  to  reproduce  itself,  it  is  necessary 


that 


r-' 


(3-24) 


Comparison  of  (3*23)  and  (3*24)  with  (3-l6)  and  (3*17)  shows  that 
formulating  the  problem  in  terms  of  the  finite  difference  approxi- 
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mation  to  the  adjoint  equation  (3*1*0  leads  to  the  same  set  of 
equations  as  the  method  of  Nelson  and  De  Marcus. 


IV.  MODERATION  BY  HYDROGEN 


In  this  section  the  general  treatment  of  Section  II  will  be 

applied  to  a  reactor  containing  hydrogen  and  an  element  which  scatters 

isotropically  and  without  energy  loss^«  It  is  convenient  to  use  the 

neutron  energy  as  an  independent  variable  instead  of  the  velocity, 

since  the  distribution  in  energy  of  neutrons  scattered  by  hydrogen 

2 

is  uniform  and  depends  only  on  the  initial  energy  . 

If  n(E)dE  is  re-defined  as  the  rate  at  which  neutrons  leave 
the  energy  interval  dE  by  any  means,  and  similar  definitions  are 
made  for  the  other  quantities,  equations  (2.17)  and  (2.18)  become 


m(e) 


.  [  de‘  /**(£')  +■  J<A£'Ka<?eO /»<*')  (4#1) 

J  c  r 


where  E0  is  the  energy  corresponding  to  v.  Since  =  2/3  for 
•5 

hydrogen  ,  the  kernels  are 


Similar  results  have  "been  obtained  by  J.  Ashkin  at  Los  Alamos 
under  these  assumptions,  and  with  a  more  general  scattering  law 
for  hydrogen.  I  am  Indebted  to  T.  A.  Welton  for  bringing  this 
to  my  attention. 

2  H.  Soodak  and  E.  C.  Campbell,  Elementary  Pile  Theory,  AECD-2201, 
Oak  Ridge,  Technical  Information  Branch,  Atomic  Energy  Commission, 
May  19^9,  P*  5* 
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^'Ve)  - 

cr 

\ 

e 

a  1^. 

cA 

°s« 

<r 

t 

e 

l-fi 

a~sn 

» 

ou 

<r 

e 

= 

ft  H  J_ 

~  c 

where  csh  is  the  scattering  cross-section  of  the  hydrogen  present, 
since  only  collisions  with  hydrogen  result  in  loss  of  energy. 
Differentiating  (4.1)  and  (4.2)  with  respect  to  E  yields  the 
following  simultaneous  differential  equations: 

As 

(4.4) 

^vw  tE)  -  klc4Ve)  /*(€) 

<U  '  ** 

provided  X(E)  is  differentiable,  with  the  boundary  conditions 


(et)  ~  0 

l%)  -  X  teJ 


(4.5) 


obtained  from  (It- .1)  and  (It. 2)  by  letting  E  approach  EQ. 

Instead  of  using  this  general  source  term,  however,  we  will 
solve  (4.1)  and  (4.2)  for  the  case 

X  <E)  =  S  lE-G o)  (1^6) 

where  6  (E-E0)  is  the  Dirac  delta  function,  corresponding  to  mono- 
energetic  fission  neutrons. 
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Let 

-  a!  (E )  +  <^Ce-e0>  (4*7) 

where  “n(E)  is  a  bounded  function*  Substituting  this  together  with 
(4.6)  in  (4.1)  and  (4.2), 

^  (6^  «  j  die'  K0>(€»)  /wi  IE')  J^-e‘  fca\e*)  &  (SO  + 

-  (E)  =  -  j^le'  kaV>  w(£‘)  +  L*  &(*')%  (e‘)  (U,8) 

E  **  -Hc^eo 

Letting  E  approach  EQ,  the  boundary  conditions  on  m  and  n  are 

0)  -  KrtVee) 

x  ^  )  (4-9) 

/h  (e0)  *  ^ 

In  order  to  obtain  an  explicit  solution  we  assume  constant 
cross-section,  so  that  using  (4.3)  the  kernels  can  be  written 

K0)(E)  -  (4.10) 


Differentiating  (4.8)  leads  to 


de 

*  f!  /m  te) 

£ 

Si 

C 

<Ov  (£) 

«4vC 

de 

£l  /v*t  (€) 
E 

- 

£ 

AV  (€) 

We  assume  solutions  of  the  form 


/y* t-  (E) 

^  (€) 


(4.11) 


(4.12) 
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where  the  constants  a,  b  and  p  are  to  be  determined.  Substituting 
in  (^.11)  leads  to  the  equation 

)  SV**  +  c*  -*=T  o 


-  c3  +■  (f>  +£$)'*  -  o 


(4.13) 


(4.18) 
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Solving  ( il- . IT )  and.  (4.18)  for  aA >a2,  bj  and  bg,  the  solution  to 
(4.1)  and  (4.2)  is 


An  LE)  =  - 


c ,  cx  4C-*  C4  +  Ci  Pj, 


Cp,  -  PO  e, 


(£)• 


4-  Cj  Pj 

Cfi  -p>  )£c 


^llE)  =  C,C7  -C„C4  -dPy  /£  y  ‘  -C3C^4.C4C44C,Pt 

cp,  -  o  ec  *  (p,  -  PO  e0 


Substituting  these  1&  (2. 30)  and  (2*31)  yields  expressions  for 


the  flux  and  current- 
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